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I describe a subtraction scheme for the next-to-next-to-leading order calculation of single inclusive production 
at hadron colliders. Such processes include Drell-Yan, W^, Z and Higgs Boson production. The key to such 
a calculation is a treatment of initial state radiation which preserves the production characteristics, such as the 
rapidity distribution, of the process involved. The method builds upon the Dipole Formalism and, with proper 
modifications, could be applied to deep inelastic scattering and e+e^ annihilation to hadrons. 
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INTRODUCTION 

Next-to-next-to-leading order (NNLO) calculations com- 
bine three contributions: Second order virtual corrections, 
first order virtual corrections to single real emission and dou- 
ble real emission. For special processes, like inclusive sin- 
gle particle production or inclusive deep inelastic scattering, 
fully inclusive NNLO calculations can be performed by ana- 
lytically integrating over the total phase space [1-9]. In doing 
so, however, much of the exclusive information from the real 
emission processes is lost. In addition, it is difficult to impose 
geometric and kinematic cuts on the phase space that would 
correspond to a realistic experimental environment, although 
important progress within the framework of totally inclusive 
calculations has recently been made on this subject [10, 11]. 
Still, one would like to be able to perform a numerical calcu- 
lation, valid to NNLO, which would retain the exclusive infor- 
mation of hard real emission and allow one to impose arbitrar- 
ily complicated cuts on the data. 

The difficulty in constructing such a program is that each of 
the three contributions to the NNLO cross section is infrared 
divergent. It is only the sum of the three that yields a mean- 
ingful physical result. Combining the three terms is made dif- 
ficult by the fact that each term involves a different number 
of final state particles and must therefore be integrated over 
a different phase space. Still, it should be possible to reorga- 
nize the infrared singularities, by adding and subtracting local 
counter-terms, so that one can perform and combine three ^- 
nite calculations. Indeed, at next-to-leading order (NLO), the 
solution is well understood and a variety of methods have been 
successfully implemented. 

A great deal of work has been devoted to the develop- 
ment of a general algorithm for constructing the local counter- 
terms at NNLO [12-17], but as yet there is no concrete ex- 
ample of a successful algorithm. Indeed the first successful 
"semi-inclusive" NNLO calculation was accomplished very 
recently [18], and did not use the analytic cancellation of in- 



frared singularities afforded by a general subtraction method. 
Instead, it used the strategy of sector decomposition [19-23] 
to numerically cancel the singular contributions. 

In this letter, I will describe an analytic subtraction scheme 
for the NNLO calculation of inclusive single particle produc- 
tion which builds upon the framework of the dipole formal- 
ism [24, 25]. In doing so, I do not solve the hardest problems 
associated with constructing a general purpose NNLO algo- 
rithm. Instead, I make use of the relative simplicity of the 
amplitudes involved and the control afforded by my ability to 
compute total integrals in each of the phase spaces. 

THE BASIC FRAMEWORK AT NLO 

I first review the basic framework for NLO calculations and 
build upon that to construct a solution at NNLO. An NLO cal- 
culation for an n parton cross section can be expressed as 

a^^O^f daf},+ [da^^\ (1) 

where the subscripts refer to the number of final state particles 
for a particular term and the superscripts refer to the order in 
the a.s expansion to which the term is calculated. Because 
of the different numbers of particle in their phase spaces, the 
two calculations must be performed separately. The fact that 
both contributions are infrared divergent presents a non-trivial 
challenge to numerical calculations. A solution to this prob- 
lem involves the construction of local counter-terms to the 
real emission cross section which permit a factorization of the 
phase space. The construction of such local counter-terms is 
made possible by the infrared factorization properties of QCD 
matrix elements. Introducing such terms, the NLO calculation 
now becomes 

Jn 7/1-1-1 
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where daj^J^ is the local counter-term. daj^^J^ is constructed 

such that it factorizes into the product of t/cri"', the n parton 
Born term, and an infrared singular term that can be com- 
pletely integrated out of a factorization of (n + l)-body phase 
space into n-body phase space times 1-body phase space. 
Thus, the terms on the first line in Eq. (2) are evaluated in 
the (n + l)-body phase space, while those on the second line 
are evaluated in the n-body phase space. 

When the radiative emission is hard, or rather, resolved, 
daj^^^ and dajl^^ map to different points in phase space and 
may be quite different in numerical value. When the emission 
is unresolved, however, the two terms map to the same point 
in phase space and their infrared singular terms cancel numer- 
ically. In the limit, the cancellation is complete and such terms 
do not contribute to the integral. 
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In analogous fashion, an NNLO calculation is expressed as 



^NNLO 



n+2 Jn+L 



doP. (3) 



Again, each of the integrals on the right hand side is in- 
frared divergent and again, one proceeds by constructing local 
counter-terms, but now the structure of the terms will be rather 
more complicated. 

aNNLO^I 

Jn+2 ^ ' 

/ + / (.afi-^e) (4) 
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where daj^^2 is the tree-level single real emission counter- 
term to dGj^^2 lil^s that found at NLO, dpf^2 ^ local counter- 
term to da,^^2 which cancels the singularities due to a sub- 
sequent second real emission, dy^^2 the tree-level double 
real emission counter-term to daf^2' ^'^i+i 
loop single real emission counter-term to ddj^^^-^ . As I have 
written Eq. (4), the terms on the first line are all computed in 
the {n + 2)-body phase space, the terms on the second line in 
(n -|- l)-body phase space (after integrating out a single emis- 
sion from daj^2 d^^^j) the terms on the third hne 
are computed in n-body phase after integrating out the single 
and double emissions from da^^^ and d'Y^^2 respectively. 

It might seem that dpj-^2 could be constructed by iterat- 
ing the procedure used to produce da^^^2- This is unlikely to 
work properly since only the infrared structure of each fac- 
torization is universal. The overlap of infrared singularities 
from one factorization with the finite remainder of the other 
is likely to generate spurious infrared divergences. Moreover 



it is not necessary. One does not need a factorized approxi- 
mation to da^2 since one does not need to integrate out the 
second emission to get down to the n-body phase space. In- 
stead, one merely needs to know the locus in phase space that 



daj^2 will map to under sequential emission. Thus, dj3^ 
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should be made numerically identical to daj^2 t'ut should be 
evaluated at the point in phase space corresponding to the sec- 
ond emission. 

The dYjj^2 ^^i+i counter-terms must approximate the 
soft and collinear limits of one-loop single real emission [26- 
28] and double real emission [29-33]. Recently, Weinzierl has 

reported a result for t/o:^^\ appropriate for one-loop final-state 
emission, which is the only case needed for the computation 
of e+e" — > jets [16]. 

It would be possible to reformulate Weinzierl's subtraction 
terms for case of initial state radiation required for the single 
particle production processes considered here, but that is not 
the solution that I propose. Instead, I use the same strategy 
in constructing dYj^j2 ™d <io:^^'j as 1 used for constructing 

rather than constructing approximations to the matrix 
element at some point in phase space, I use the exact matrix 
elements but evaluate them as if they were at different point 
in phase space. I can do this because 1 have complete an- 
alytic control of the total integrals of these matrix elements 
over phase space. 

With this strategy, it is clear that if I compute the total rate, 
placing no kinematic or geometric cuts on the configuration 
and not binning any distributions, the first two lines in Eq. (4) 
vanish identically, while the third hne gives the known result 
for the inclusive cross section. 



RAPIDITY DISTRIBUTIONS 

The simple framework described above is sufficient for de- 
scribing the total rate, but does not include all of the informa- 
tion available in inclusive production. In addition to the total 
rate, one can also observe the rapidity distributions of the vec- 
tor or Higgs boson [10, 11]. In order to reproduce rapidity 
distributions in this calculation, one needs rapidity informa- 
tion in the subtraction term. The subtraction term in the full 
(n + 2)-body phase space certainly contains this information. 
What is needed is a means of capturing this information in the 
n-body phase space. The way to do so is to understand the 
structure of the subtraction term in the «-body phase space: in 
final state emission, the subtraction term in the n-body phase 
space is essentially a number; in initial state emission, it is 
a convolution of the emission contributions with the partonic 
cross section and can be expressed in a form reminiscent of 
mass factorization. 
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where the superscripts refer to the order in the expansion in 
Us of each term. 

Unlike mass factorization, which relates the (infrared sin- 
gular) contributions determined by sum of squared matrix el- 
ements integrated over phase space to a convolution of finite 
partonic cross section with the mass factorization counter- 
terms associated with the Altarelli-Parisi splitting functions, 
Eq. (5) seeks to relate individual components, say qq Vgg, 
to a convolution of the virtual cross section, evaluated to the 
appropriate order, with real emission counter- terms. The "ap- 
propriate order" is determined by the factorization properties 
of the QCD matrix elements. Specifically, double real ra- 
diation terms map onto a convolution of the virtual term at 
Born level with either two single emission terms or one dou- 
ble emission term, while one-loop single real radiation terms 
map onto the sum of the one-loop virtual term convolved with 
one first-order single emission term and the Born level virtual 
term convolved with one second-order single emission term. 

To be more explicit, 1 will consider vector boson produc- 
tion and define the mass of the vector boson to be My and the 
center-of-mass (cm.) energy squared of the production pro- 
cess to be S. The ratio of the vector boson mass squared to 
s is defined to be z = My / s, the momentum fractions of the 
incoming partons are defined to be xi and X2 (s = x\ X2 s), and 
the fraction of those momentum fractions that go into vector 
boson production are defined to be wi and W2 (z = wi wi). In 
terms of these parameters, 

(cc / cc \ ^ \ 
ao + — ai + (— a2 , 
n \n) ) 

fKwO = 5(1 - wO + |fW(wO + (f )'ff . 

In the virtual terms, all of the energy goes into vector boson 
production and the vector boson rapidity is identical to the 
rapidity of the lab system. 
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In the case of single real radiation, some of the cm. energy 
goes into vector boson production, but some goes into a sin- 
gle collinear emission. The convolution consists of terms for 
which the emission is in the direction of parton 1 and terms 
for which it is in the direction of parton 2. Given the energy 
and boost of the cm. system, the mass of the vector boson and 
the direction of the radiative emission, it is a simple exercise 
to determine the rapidity of the vector boson: 



(8) 



In an actual calculation, one generates vector bosons at a given 
rapidity and then convolves the parton distributions over the 
parton momentum fractions, using f to properly weight the 
contributions above threshold. 

Things are more complicated when one must consider dou- 
ble emission. In the triple collinear Umit, where both emis- 
sions come from the same incoming leg, the kinematics is un- 
changed from the single emission case and one can determine 



the rapidity of the vector boson from the direction of the emis- 
sion. The case of overlapping emission, where each incoming 
parton contributes to the radiation is more complicated. In or- 
der to determine the rapidity of the vector boson, one needs 
to know what fraction of the radiation came from each side. 
Thus, for given parton momentum fractions xi , X2, rather than 
generating vector bosons at distinct rapidities, as in Eq. (8), 
one generates them over a continuum of rapidities: 



1 XlWl 



2 X2,W2 

-In— -|--lnz < yn-^i < ^In ^Inz. 

2 X2 2 2 X2 2 



(9) 



Again, one generates vector bosons of a given rapidity and 
then convolves the parton distributions over the momentum 
fractions using F to properly weight the contributions above 
threshold, but now one needs to know the appropriate F, on 
each side. The total integral over double emission phase space 
does not give f however, it gives the second order contribu- 
tion to the convolution 

f = rf\z) + + f ^')(wi) ® f «(W2) . (10) 

~ (2) 

The terms correspond to the triple collinear limits of 
the incoming partons. Subtracting these terms off, one is left 
with just the overlapping emission term, each component of 
which can be expanded as a Laurent series in e, 

overlap " ^ 1 <^ 2 ' 



p{2) ^ £(-4) 8(-3l £(-2) g(-l) 
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Since the infrared structure of single emission is universal, the 
most singular terms in Eq. (11) are fixed and the less singular 
terms can be solved for, term by term. When single emission 
is considered in isolation, the nonsingular terms, aj;„>o, are 
non-universal. In the context of overlapping divergences from 
double emission, however, these terms can be given meaning- 
ful definitions. 

g(-4) = «l,(-2) a2,(-2) I 

8{-3) = fll,(-l) «2,(-2) +flL(-2) a2,(-l) > 
^(-2) = ai,(0) ® «2,(-2) + ai,(-2) <8) 02,(0) 

+ ai,(-i)® a2,(-i)> 
g(_i) = ai_(i) (g) a2,(_2) + ai,(_2) a2,(i) (12) 

+ «1,(0) ® «2,(-l) <8> 02,(0) ) 

^(0) = «1,(2) ^ «2,(-2) +«l,(-2) O «2,(2) 

+ «1,(1) ® a2,(-l) +«1,(-1) <8> «2,(-l) 
+ «1,(0)® «2,(0)> 

The singular terms must cancel against other infrared contri- 
butions, so all components involved in those terms can be de- 
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termined unambiguously. While the soft contributions (cor- 
responding to delta functions in the parton fractions) cannot 
be uniquely distributed among the separate terms, their distri- 
bution does not affect the magnitude of the total subtraction 
nor its location in phase space. The only room for ambiguity 
would come in the terms that are new to ^(o), i.e. «i,(2) and 
02.(2) ■ When the splittings are symmetric, say for qq Vgg 
each splitting would bcq^ qg (where d identifies the incom- 
ing parton), the solution is symmetric and there is no ambi- 
guity. When the splittings are asymmetric, say for qg Vqg 
where one splitting is q ^ qg and the other is q ^ gq (note 
the order), it may seem that one cannot separate the new terms. 
However, the prior knowledge of the more singular terms with 
which they are convolved should allow for this extraction. 
Certainly the dominant terms can be determined unambigu- 
ously. It is also possible to obtain rapidity distributions at 
NNLO within the framework of inclusive calculations and this 
has already been done for Drell-Yan production [11]. The in- 
formation from such a calculation could be used to resolve 
any lingering ambiguities. 

Imphcit to the preceding discussion is that the ^^'^ 

<ia^^\ in the (n -|- 2) and (n -|- l)-body phase spaces must be 
evaluated at points in phase space that correspond to those de- 
scribed in n-body phase space. Those points correspond to 
conserving the incoming parton momentum fractions of the 
true process, but mapping all final state momenta onto the 
beam axes. In the event of unresolved emission, this mapping 
ensures a proper subtraction of the singularities. 



CONCLUSIONS 

In this letter, I have described the explicit construction of 
the subtraction terms needed for an NNLO parton-level monte 
carlo calculation of single inclusive production (Higgs bosons, 
vector bosons, Drell-Yan) within the dipole formalism. For 
these simple processes, one has complete analytic control over 
the total integrals of the matrix elements. This control allows 
one to use the matrix elements themselves as the subtraction 
terms. In order to capture the full information available from 
inclusive calculations, one needs to retain information on both 
the total rate and the rapidity distribution of production. I have 
outlined a method for doing so. 
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Note: As this paper was being completed, I became aware 
of Ref. [34], in which essentially the same solution is pro- 
posed for the final state emission process, e~^e~ — > 2 jets, and 
is worked out in great detail. 
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